#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "fft_math.h"

#define PI 3.1415926535
#define N 1024 
#define M 10
	
void printBin(unsigned i, int bits) {
    unsigned  shiftcount = bits;
    while(bits-- != 0)
        i & ( 1 << bits ) ? printf("1") : printf("0");
}
void main(void)
{
	int i,j,k,r;
	int p,L,B;
	unsigned int I,J,K,F0,F1,m,n;
	float SinIn[N];
	float dataR[N],dataI[N];
	float dataA[N];
	float Tr,Ti,temp;
	FILE* fp;
	int f1 = 100;//   	# 信号的特征频率1
	int f2 = 300;//   # 信号的特征频率2
	int f3 = 500;//   # 信号的特征频率3
	printf("[%s +%d %s] n=%d log2(N)=%f\n", __FILE__, __LINE__, __func__, N, log2(N));
	printf("[%s +%d %s] f1=%d f2=%d f3=%d\n", __FILE__, __LINE__, __func__, f1, f2, f3);

	for(i = 0; i < N; i++)
	{
		float j = i;
		SinIn[i] = 10 + 7 * sin(2*PI*f1*(j*1/N)) + 5 * sin(2*PI*f2* (j*1/N)) + 3 * sin(2*PI*f3* (j*1/N) );//A0 * (sin(2*PI*i/25)+sin(2*PI* i * 0.4 ));
		dataR[i] = SinIn[i];
		dataI[i] = 0;
		dataA[i] = 0;
	}
	if((fp=fopen("input.txt","wb"))==NULL)
	{
 		printf("fopen input.txt failed!\n");
 		exit(0);
	}
	for(i = 0; i < N ; i++)
		fprintf(fp,"%f,",dataR[i]);
	fclose(fp);
	FFT(dataR,dataI,dataA,N,M);
	for(i=0;i<N;i++)
		printf("%f,",dataA[i]);
	if((fp=fopen("out_fft.txt","wb"))==NULL)
	{
 		printf("fopen out_fft.txt failed!\n");
 		exit(0);
	}
	for(i = 0; i < N ; i++)
		fprintf(fp,"%f,",dataA[i]);
	//fwrite(&dataR,sizeof(dataR),1,fp);
	fclose(fp);
	for(i = 0; i < N; i++)
	{
		float j = i;
		SinIn[i] = 10 + 7 * sin(2*PI*f1*(j*1/N)) + 5 * sin(2*PI*f2* (j*1/N)) + 3 * sin(2*PI*f3* (j*1/N) );//A0 * (sin(2*PI*i/25)+sin(2*PI* i * 0.4 ));
		dataR[i] = SinIn[i];
		dataI[i] = 0;
		dataA[i] = 0;
	}
	FFTR(dataR,dataI,dataA,N,M);
	for(i=0;i<N;i++)
		printf("%f,",dataA[i]);
	if((fp=fopen("out_fftr.txt","wb"))==NULL)
	{
 		printf("fopen out_fftr.txt failed!\n");
 		exit(0);
	}
	for(i = 0; i < N ; i++)
		fprintf(fp,"%f,",dataA[i]);
	//fwrite(&dataR,sizeof(dataR),1,fp);
	fclose(fp);
	for(i = 0; i < N; i++)
	{
		float j = i;
		SinIn[i] = 10 + 7 * sin(2*PI*f1*(j*1/N)) + 5 * sin(2*PI*f2* (j*1/N)) + 3 * sin(2*PI*f3* (j*1/N) );//A0 * (sin(2*PI*i/25)+sin(2*PI* i * 0.4 ));
		dataR[i] = SinIn[i];
		dataI[i] = 0;
		dataA[i] = 0;
	}
	printf("[%s +%d %s] N=%d\n", __FILE__, __LINE__, __func__, N);
	printf("[%s +%d %s] M=%d\n", __FILE__, __LINE__, __func__, M);
	for(I=0;I< N;I++)
	{  
		J=0;
		for(k=0;k<(M/2+0.5);k++)
		{
			m=1;
			n=(unsigned int)pow(2,M-1);
			//printf("[%s +%d %s] M=%d k=%d\n", __FILE__, __LINE__, __func__, M, k);

			m <<= k;
			n >>= k;
			F0=I & n;
			F1=I & m;
			if(F0) J=J | m;
			if(F1) J=J | n;
		}

		if(I<J)
		{
			temp = dataR[I];
			dataR[I] = dataR[J];
			dataR[J] = temp;
		}                                
	} 


	for(L=1; L<=M;L++)
	{		
		B = 1;
		B = (int)pow(2,L-1);	
		//printf("[%s +%d %s] M=%d L=%d\n", __FILE__, __LINE__, __func__, M, L);
		for(j=0;j<=B-1;j++)		
		{
			k=1;
			k = (int)pow(2,M-L);
			//printf("[%s +%d %s] M=%d j=%d\n", __FILE__, __LINE__, __func__, M, j);
			p=1;
			p=j*k;
			for(i=0; i<=k-1;i++) 
				{
					printf("[%s +%d %s] M=%d i=%d\n", __FILE__, __LINE__, __func__, M, i);
					r=1;
					r=j+2*B*i;
					Tr=dataR[r+B]*cos(2.0*PI*p/N) + dataI[r+B]*sin(2.0*PI*p/N);
					Ti=dataI[r+B]*cos(2.0*PI*p/N) - dataR[r+B]*sin(2.0*PI*p/N);
					dataR[r+B]=dataR[r]-Tr;
					dataI[r+B]=dataI[r]-Ti;
					dataR[r]=dataR[r]+Tr;
					dataI[r]=dataI[r]+Ti;
				}	
			}
	}

	for ( i=0;i<N;i++ )
	{ 		
	 	dataA[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
	}
	if((fp=fopen("out_main.txt","wb"))==NULL)
	{
 		printf("fopen out_main.txt failed!\n");
 		exit(0);
	}
	for(i = 0; i < N ; i++)
		fprintf(fp,"%f,",dataA[i]);
	//printf("%f\n",dataA[i]);}

	//fwrite(&dataR,sizeof(dataR),1,fp);
	fclose(fp);
	printf("\n");
}

